網頁應用程式4

林嶔 (Lin, Chin)

Lesson 15

第一節 Google圖表在網頁應用程式的應用(1)

– 學員可以在googleVis examples找到一些例子

– 請複製貼上以下範例在Rstudio中

– 請至這裡下載範例資料

dat = read.csv("Example data.csv", header = TRUE)
head(dat)
##       eGFR Disease Survival.time Death Diabetes Cancer      SBP      DBP
## 1 34.65379       1     0.4771037     0        0      1 121.2353 121.3079
## 2 37.21183       1     3.0704424     0        1      1 122.2000 122.6283
## 3 32.60074       1     0.2607117     1        0      0 118.9136 121.7621
## 4 29.68481       1            NA    NA        0      0 118.2212 112.7043
## 5 28.35726       0     0.1681673     1        0      0 116.7469 115.7705
## 6 33.95012       1     1.2238556     0        0      0 119.9936 116.3872
##   Education Income
## 1         2      0
## 2         2      0
## 3         0      0
## 4         1      0
## 5         0      0
## 6         1      0
library(googleVis)
TAB = table(dat$Income)
TAB = data.frame(TAB)
colnames(TAB) = c("Income", "Freq")
TAB[,1] = c("Low income", "Middle-class", "Wealthy")
Pie = gvisPieChart(TAB)
plot(Pie)

第一節 Google圖表在網頁應用程式的應用(2)

F3_1

library(survival)
data(ovarian)    #呼叫ovarian dataset
dat = ovarian    #將ovarian轉存為dat
model <- coxph(Surv(futime, fustat) ~ age, data = dat) #利用age預測存活時間
summary(model)   #看結果
## Call:
## coxph(formula = Surv(futime, fustat) ~ age, data = dat)
## 
##   n= 26, number of events= 12 
## 
##        coef exp(coef) se(coef)     z Pr(>|z|)   
## age 0.16162   1.17541  0.04974 3.249  0.00116 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##     exp(coef) exp(-coef) lower .95 upper .95
## age     1.175     0.8508     1.066     1.296
## 
## Concordance= 0.784  (se = 0.091 )
## Rsquare= 0.423   (max possible= 0.932 )
## Likelihood ratio test= 14.29  on 1 df,   p=0.0001564
## Wald test            = 10.56  on 1 df,   p=0.001157
## Score (logrank) test = 12.26  on 1 df,   p=0.0004629

第一節 Google圖表在網頁應用程式的應用(3)

– 因此我們可以透過前後相減，獲得每個時間點的所增加的hazard

h0.hazard = basehaz(model,centered=FALSE) #取得baseline累積hazard
h0.hazard$hazard.lag = c(0,h0.hazard$hazard[-nrow(h0.hazard)]) #產生前一次的值
h0.hazard$hazard.dif = h0.hazard$hazard-h0.hazard$hazard.lag #與前一次相減
h0.hazard = h0.hazard[,c(2,4)]   #僅留下time & hazard.dif
head(h0.hazard)
##   time   hazard.dif
## 1   59 1.376452e-06
## 2  115 1.647152e-06
## 3  156 2.284238e-06
## 4  268 2.554178e-06
## 5  329 4.506719e-06
## 6  353 4.528477e-06
# 為了在Google chart中呈現每個時間點的累積機率，我們必須執行以下動作
h0.hazard.extend = data.frame(time = 0:max(h0.hazard[,1]), hazard.dif = 0) 
h0.hazard = h0.hazard[h0.hazard[,2]!=0,]
h0.hazard = rbind(h0.hazard.extend,h0.hazard)
h0.hazard = h0.hazard[order(h0.hazard[,1]),] #將每個時間點都納入Table中
head(h0.hazard)
##   time hazard.dif
## 1    0          0
## 2    1          0
## 3    2          0
## 4    3          0
## 5    4          0
## 6    5          0

第一節 Google圖表在網頁應用程式的應用(4)

NEW <- as.matrix(data.frame(age=50)) #產生一個新的資料表描述自訂參數，並儲存成矩陣格式
matrix.coef = matrix(model$coef,nrow=length(model$coef),ncol=1) #將Cox model的係數轉存為一個矩陣
indv.hazardratio = exp(NEW%*%matrix.coef) #取得個體的hazard ratio(%*%為矩陣乘法)
indv.hazard = indv.hazardratio*h0.hazard$hazard.dif #取得個體在每一個時間點的hazard
indv.cumhazard = cumsum(indv.hazard) #取得個體在各時間點的『累積』hazard
indv.cumrate = exp(-indv.cumhazard)  #轉換為累積存活率
Predic.Survival = data.frame(time = h0.hazard$time, rate = indv.cumrate) #產生預測的時間&累積存活率資料表
Predic.Survival[,2] = round(Predic.Survival[,2]*100,2) #將累積存活率轉換成百分率
head(Predic.Survival)
##   time rate
## 1    0  100
## 2    1  100
## 3    2  100
## 4    3  100
## 5    4  100
## 6    5  100

第一節 Google圖表在網頁應用程式的應用(5)

– 註：使用者可以利用help(gvisScatterChart)微調圖形參數

library(googleVis)
Scatter <- gvisScatterChart(Predic.Survival, 
                            options=list(
                              explorer="{actions: ['dragToZoom', 
                                        'rightClickToReset'],
                                        maxZoomIn:0.05}",
                              legend="none",
                              lineWidth=2, pointSize=0,
                              vAxis="{title:'Survival (%)'}",
                              vAxes="[{viewWindowMode:'explicit',
                                     viewWindow:{min:0, max:100}}]",
                              hAxis="{title:'Time (days)'}", 
                              colors="['#ff0000']",
                              width=800, height=700))
plot(Scatter)

第一節 Google圖表在網頁應用程式的應用(6)

– ui.R

library(shiny)
library(survival)
library(googleVis)

fluidPage(
  sliderInput("Age", "Please enter your age", min=40, max=80, value=50),
  htmlOutput("chart1")
)

– server.R

library(shiny)
library(survival)
library(googleVis)

######################################################
# 這些函數只需要跑1次即可

data(ovarian)   
dat = ovarian 
model <- coxph(Surv(futime, fustat) ~ age, data = dat) 
h0.hazard = basehaz(model,centered=FALSE) 
h0.hazard$hazard.lag = c(0,h0.hazard$hazard[-nrow(h0.hazard)])
h0.hazard$hazard.dif = h0.hazard$hazard-h0.hazard$hazard.lag
h0.hazard = h0.hazard[,c(2,4)]
h0.hazard.extend = data.frame(time = 0:max(h0.hazard[,1]), hazard.dif = 0) 
h0.hazard = h0.hazard[h0.hazard[,2]!=0,]
h0.hazard = rbind(h0.hazard.extend,h0.hazard)
h0.hazard = h0.hazard[order(h0.hazard[,1]),]
matrix.coef = matrix(model$coef,nrow=length(model$coef),ncol=1)

######################################################

shinyServer(function(input, output, session) {
  
  output$chart1<- renderGvis({
    NEW <- as.matrix(data.frame(age=input$Age)) #年齡=input$Age
    
    indv.hazardratio = exp(NEW%*%matrix.coef) 
    indv.hazard = indv.hazardratio*h0.hazard$hazard.dif
    indv.cumhazard = cumsum(indv.hazard)
    indv.cumrate = exp(-indv.cumhazard)
    Predic.Survival = data.frame(time = h0.hazard$time, rate = indv.cumrate) 
    Predic.Survival[,2] = round(Predic.Survival[,2]*100,2)
    
    Scatter <- gvisScatterChart(Predic.Survival, 
                                options=list(
                                  explorer="{actions: ['dragToZoom', 
                                  'rightClickToReset'],
                                  maxZoomIn:0.05}",
                                  legend="none",
                                  lineWidth=2, pointSize=0,
                                  vAxis="{title:'Survival (%)'}",
                                  vAxes="[{viewWindowMode:'explicit',
                                  viewWindow:{min:0, max:100}}]",
                                  hAxis="{title:'Time (days)'}", 
                                  colors="['#ff0000']",
                                  width=800, height=500))
    Scatter
  })
  
})

練習-1

– 請將剛剛的WebApp改寫，讓使用者能輸出age+rx的數值並進行預測

– 注意，rx的數值僅可以是1或2，請用radioButtons()來讓使用者輸入參數

– 註：radioButtons()回傳的物件為『文字』，需使用as.numeric()來使該物件轉換為數字

data(ovarian)   
dat = ovarian 
model <- coxph(Surv(futime, fustat) ~ age + rx, data = dat) 
h0.hazard = basehaz(model,centered=FALSE) 
h0.hazard$hazard.lag = c(0,h0.hazard$hazard[-nrow(h0.hazard)])
h0.hazard$hazard.dif = h0.hazard$hazard-h0.hazard$hazard.lag
h0.hazard = h0.hazard[,c(2,4)]
h0.hazard.extend = data.frame(time = 0:max(h0.hazard[,1]), hazard.dif = 0) 
h0.hazard = h0.hazard[h0.hazard[,2]!=0,]
h0.hazard = rbind(h0.hazard.extend,h0.hazard)
h0.hazard = h0.hazard[order(h0.hazard[,1]),]
matrix.coef = matrix(model$coef,nrow=length(model$coef),ncol=1)

NEW <- as.matrix(data.frame(age=50,rx=1)) #年齡=50,rx=1
    
indv.hazardratio = exp(NEW%*%matrix.coef) 
indv.hazard = indv.hazardratio*h0.hazard$hazard.dif
indv.cumhazard = cumsum(indv.hazard)
indv.cumrate = exp(-indv.cumhazard)
Predic.Survival = data.frame(time = h0.hazard$time, rate = indv.cumrate) 
Predic.Survival[,2] = round(Predic.Survival[,2]*100,2)

Scatter <- gvisScatterChart(Predic.Survival, 
                                options=list(
                                  explorer="{actions: ['dragToZoom', 
                                  'rightClickToReset'],
                                  maxZoomIn:0.05}",
                                  legend="none",
                                  lineWidth=2, pointSize=0,
                                  vAxis="{title:'Survival (%)'}",
                                  vAxes="[{viewWindowMode:'explicit',
                                  viewWindow:{min:0, max:100}}]",
                                  hAxis="{title:'Time (days)'}", 
                                  colors="['#ff0000']",
                                  width=800, height=500))
plot(Scatter)

練習-1 答案

library(shiny)
library(survival)
library(googleVis)

fluidPage(
  sliderInput("Age", "Please enter your age", min=40, max=80, value=50),
  radioButtons("rx", "Please select a treatment group", c("1","2")),
  htmlOutput("chart1")
)
library(shiny)
library(survival)
library(googleVis)

######################################################
# 這些函數只需要跑1次即可

data(ovarian)   
dat = ovarian 
model <- coxph(Surv(futime, fustat) ~ age + rx, data = dat) 
h0.hazard = basehaz(model,centered=FALSE) 
h0.hazard$hazard.lag = c(0,h0.hazard$hazard[-nrow(h0.hazard)])
h0.hazard$hazard.dif = h0.hazard$hazard-h0.hazard$hazard.lag
h0.hazard = h0.hazard[,c(2,4)]
h0.hazard.extend = data.frame(time = 0:max(h0.hazard[,1]), hazard.dif = 0) 
h0.hazard = h0.hazard[h0.hazard[,2]!=0,]
h0.hazard = rbind(h0.hazard.extend,h0.hazard)
h0.hazard = h0.hazard[order(h0.hazard[,1]),]
matrix.coef = matrix(model$coef,nrow=length(model$coef),ncol=1)

######################################################

shinyServer(function(input, output, session) {
  
  output$chart1<- renderGvis({
    NEW <- as.matrix(data.frame(age=input$Age,rx=as.numeric(input$rx))) #年齡=input$Age;組別=as.numeric(input$rx)
    
    indv.hazardratio = exp(NEW%*%matrix.coef) 
    indv.hazard = indv.hazardratio*h0.hazard$hazard.dif
    indv.cumhazard = cumsum(indv.hazard)
    indv.cumrate = exp(-indv.cumhazard)
    Predic.Survival = data.frame(time = h0.hazard$time, rate = indv.cumrate) 
    Predic.Survival[,2] = round(Predic.Survival[,2]*100,2)
    
    Scatter <- gvisScatterChart(Predic.Survival, 
                                options=list(
                                  explorer="{actions: ['dragToZoom', 
                                  'rightClickToReset'],
                                  maxZoomIn:0.05}",
                                  legend="none",
                                  lineWidth=2, pointSize=0,
                                  vAxis="{title:'Survival (%)'}",
                                  vAxes="[{viewWindowMode:'explicit',
                                  viewWindow:{min:0, max:100}}]",
                                  hAxis="{title:'Time (days)'}", 
                                  colors="['#ff0000']",
                                  width=800, height=500))
    Scatter
    })
  
})

第二節 在網頁應用程式中的高耗時程序(1)

– 這筆data是來自於bioconductor/golubEsets套件中的Golub_Merge data set

data = read.csv("microarraydata.csv")     #讀檔
Y = substr(colnames(data)[-1],1,3)        #把不同類型的癌症分開
X = data[,-1]                             #取得Gene expression matrix
result_t.test = rep(NA,nrow(X))  #產生一個向量儲存所有probe的p value
for (i in 1:nrow(X)) {           #利用for迴圈來做t test (大約需要15-20秒)
  result_t.test[i] = t.test(as.numeric(X[i,])~Y)$p.value
}
result_U.test = rep(NA,nrow(X))  #產生一個向量儲存所有probe的p value
for (i in 1:nrow(X)) {           #利用for迴圈來做U test (大約需要15-20秒)
  result_U.test[i] = wilcox.test(as.numeric(X[i,])~Y,exact=FALSE)$p.value
}

第二節 在網頁應用程式中的高耗時程序(2)

– reactive()、eventReactive()可以幫助我們做到這個功能

– 另外，我們再介紹進度條函數:withProgress()，這可以提示使用者目前正在運行中。

library(shiny)

shinyUI(pageWithSidebar(
  
  headerPanel("Multiple test"), 
  
  sidebarPanel(
    fileInput(inputId="files", label=h4("Upload your data file:"), multiple=FALSE, accept="text/csv"),
    helpText("Note: you only can upload the .csv file."),
    radioButtons("method", label = h4("t test or U test?"), choices = list("t test" = "t", "U test" = "u")),
    actionButton("submit",strong("Start to analyze!"),icon("list-alt"))
  ),
  
  mainPanel(
    sliderInput("n.top", label = h4("Show first x probes:"), min = 10, max = 100, value = 50),
    downloadButton("download", label = "Download Result", class = NULL),
    tableOutput("table")
  )  
  
))
library(shiny)

shinyServer(function(input, output) {
  
  DATA <- reactive({
    if (is.null(input$files)) {return()} else {
      dat <- read.csv(input$files$datapath,header=T)
      return(dat) 
    }
  })
  
  RESULT <- eventReactive(input$submit,{
    data = DATA()
    if (is.null(data)) {return()} else {
      withProgress(message = "In processing...",value=0,{
        Y = substr(colnames(data)[-1],1,3)
        X = data[,-1]
        n.probe = nrow(X)
        result = rep(NA,nrow(X))
        for (i in 1:n.probe) {
          incProgress(1/n.probe)
          if (input$method=="t") {result[i] = t.test(as.numeric(X[i,])~Y)$p.value}
          else {result[i] = wilcox.test(as.numeric(X[i,])~Y,exact=FALSE)$p.value}
        }
        return(data.frame(probe=data[,1],p.value=result))
      })
    }
  })
  
  output$table <- renderTable({
    result = RESULT()
    if (is.null(result)) {return()} else {
      result = result[order(result[,2]),]
      return(head(result,input$n.top))
    }
  })
  
  output$download <- downloadHandler(
    filename = function() {'Result.csv'},
    content = function(con) {
      result = RESULT()
      if (is.null(result)) {return()} else {
        result = result[order(result[,2]),]
        write.csv(result,con,quote=FALSE,row.names=FALSE)
      }
    }
  )
  
})

練習-2

– 注意：ctree()不允許依變項有任何missing value，必須刪除missing後才能放入ctree函數

library(party)
data(iris)

iris = iris[!is.na(iris$Species),] # remove missing value

Result = ctree(Species~.,data=iris,
               controls=ctree_control(maxdepth=2,            # maximum depth of the tree
                                      mincriterion=0.95)     # 1-significant level
               ) 
plot(Result)

  1. 能允許使用者上傳一個典型的資料表(如範例)
  2. 請使用者分類所有的變項是連續變項或是類別變項
  3. 請使用者選擇一個變項做為依變項
  4. 請使用者決定哪些變項要進入決策樹做分析
  5. 請使用者決定樹的最大深度&顯著水準
  6. 設置一個按鈕，在按下後立即進行分析
  7. 設置一個下載鈕，讓使用者能下載Tree plot

練習-2 答案

library(shiny)
library(party)

shinyUI(pageWithSidebar(
  
  headerPanel("Tree"), 
  
  sidebarPanel(
    fileInput(inputId="files", label=h4("Upload your data file:"), multiple=FALSE, accept="text/csv"),
    helpText("Note: you only can upload the .csv file."),
    uiOutput("choose_columns1"),
    uiOutput("choose_columns2"),
    uiOutput("choose_columns3"),
    sliderInput("maxdepth", label = h4("Maximum depth of the tree:"), min = 1, max = 10, value = 2),
    numericInput("sig", label = h4("Significant level:"), min = 1e-4, max = 1-1e-4, value = 0.05),
    actionButton("submit",strong("Start to analyze!"),icon("list-alt"))
  ),
  
  mainPanel(
    downloadButton("download", label = "Download Result", class = NULL),
    plotOutput("plot")
  )  
  
))
library(shiny)
library(party)

shinyServer(function(input, output) {
  
  DATA <- reactive({
    if (is.null(input$files)) {return()} else {
      dat <- read.csv(input$files$datapath,header=T)
      return(dat) 
    }
  })
  
  output$choose_columns1 <- renderUI({ 
    dat = DATA()
    if (is.null(dat)) {return()} else {
      colnames <- colnames(dat)
      return(selectInput("cats", h4("Choose categorical variables:"), choices = colnames, multiple = TRUE))
    }
  })
  
  output$choose_columns2 <- renderUI({ 
    dat = DATA()
    if (is.null(dat)) {return()} else {
      colnames <- colnames(dat)
      selectInput("Y", h4("Choose a dependence variable:"), choices = colnames)
    }
  })
  
  output$choose_columns3 <- renderUI({ 
    dat = DATA()
    if (is.null(dat)|is.null(input$Y)) {return()} else {
      colnames <- colnames(dat)
      selectInput("X", h4("Choose independence variables:"), choices = colnames[which(colnames!=input$Y)], multiple = TRUE)
    }
  })
  
  TREE <- eventReactive(input$submit,{
    dat = DATA()
    if (is.null(dat)|is.null(input$Y)|is.null(input$X)) {return()} else {
      dat.y=data.frame(dat[,input$Y])
      if (input$Y%in%input$cats) {dat.y[,1] <- factor(dat.y[,1])}
      dat.x=data.frame(dat[,input$X])
      for (i in 1:length(input$X)) {
        if (input$X[i]%in%input$cats) {dat.x[,i]=factor(dat.x[,i])}
      }
      names(dat.x)=input$X
      dat.x=data.frame(dat.x,y1=dat.y[,1])
      dat.x=dat.x[!is.na(dat.x$y1),]
      Result = ctree(y1~.,data=dat.x,
                     controls=ctree_control(maxdepth=input$maxdepth,           
                                            mincriterion=1-input$sig)   
      )
      return(Result)
    }
  })
  
  output$plot <- renderPlot({
    result = TREE()
    if (is.null(result)) {return()} else {
      return(plot(result))
    }
  })
  
  output$download <- downloadHandler(
    filename = function() {'Tree.pdf'},
    content = function(con) {
      result = TREE()
      if (is.null(result)) {return()} else {
        pdf(con)
        plot(result)
        dev.off()
      }
    }
  )
  
})

小結

  1. conditionalPanel()
  2. 讓使用者能上傳檔案
  3. renderUI及uiOutput
  4. UI的優化
  5. Datatable & Google charts
  6. 設計server的反應流程